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In recent years, numerical solutions of the equations of compressible magnetohydrody- 
namic flows have been found to contain intermediate shocks for certain kinds of problems. 
Since these results would seem to be in conflict with the classical theory of magnetohy- 
drodynamic shocks, they have stimulated attempts to reexamine various aspects of this 
theory, in particular the role of dissipation. In this paper we study the general relation- 
ship between the evolutionary conditions for discontinuous solutions of the dissipation- 
free system and the existence and uniqueness of steady dissipative shock structures for 
systems of quasilinear conservation laws with a concave entropy function. Our results 
confirm the classical theory. We also show that the appearance of intermediate shocks 
in numerical simulations can be understood in terms of the properties of the equations 
of planar magnetohydrodynamics for which some of these shocks turn out to be evolu- 
tionary. Finally, we discuss ways in which numerical schemes can be modified in order 
to avoid the appearance of intermediate shocks in simulations with such symmetry. 



1. Introduction 

It is well known that not all discontinuous solutions of hyperbolic conservation laws are 
admissible. Some of these can be excluded on physical grounds. For example, expansion 
shocks in gas dynamics must be discarded since they do not satisfy the second law of 
thermodynamics. Others can be excluded for purely mathematical reasons such as the 
fact that they do not satisfy uniqueness and existence conditions or are structurally 
unstable with respect to small perturbations of the initial data. These mathematical 
conditions are usually called evolutionary conditions. For example, intermediate shocks 
in magnetohydrodynamics (MHD) satisfy the second law but are not evolutionary. 

This subject was extensively studied between the late 1940's and early 1960's (e.g. 
Courant & Friedrichs 1948, Lax 1957, Akhiezer et al. 1959, Germain 1960, Gel'fand 
1963, Polovin 1961) and a full account can be found in numerous textbooks (e.g. Jeffrey 
& Taniuti 1964, Cabannes 1970, Somov 1994). Until recently there was general agree- 
ment that admissible shocks must both satisfy the evolutionary condition and possess a 
steady dissipative shock structure, although the relation between these conditions was 
not entirely clear. There the matter rested until time-dependent numerical solutions of 
the dissipative MHD equations showed that certain types of intermediate shocks can 
arise from smooth initial data (Wu 1987). Shortly thereafter, Brio & Wu (1988) found 
intermediate shocks in their numerical solution for a particular MHD Riemann problem. 
More recently, intermediate shocks have been also been found in two-dimensional simula- 
tions (De Sterck et al. 1998). Furthermore, Chao et al. (1993) have reported a detection 
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of an interplanetary intermediate shock in the Voyager 1 data. All this has caused some 
authors to reject the classical theory and to suggest that the evolutionary condition is not 
relevant to dissipative MHD (Wu 1987, 1988a,b, 1990; Kennel, Blandford & Wu 1990; 
Hada 1994, Myong & Roe 1997a, b) and has led to a reexamination of the whole question 
of the existence, or otherwise, of non-classical shocks (see Glimm 1988, Freistuhler & 
Liu 1993, Myong & Roe 1997a and references therein). There are, however, others who 
argue that there is nothing wrong with the classical theory (e.g. Barmin, Kulikovsky & 
Pogorelov 1996; Falle & Komissarov 1997). 

The matter clearly needs to be resolved, particularly since the existence, or otherwise 
of intermediate shocks is of crucial importance not only for fundamental MHD processes 
such as reconnection (Wu 1995), but is also relevant to many other astrophysical applica- 
tions. The purpose of this paper is to try and clear the matter up by showing that there 
is neither a real conflict between the classical shock theory and the results of numeri- 
cal calculations nor any incompatibility between ideal and dissipative MHD. In order to 
make the discussion complete, we have put together and extended a number of results 
from the literature that have tended be ignored or misunderstood. 

This paper is organised as follows. In §|| we briefly review the classical shock theory 
and the evolutionary conditions. In §|| we study the relationship between these conditions 
and the uniqueness and existence of steady dissipative shock structures for systems with 
a concave entropy function. In §[| we apply these results to the full system of MHD 
equations and to the reduced system of planar MHD. In §[| we present the results of 
numerical calculations which show that, for both these systems, the behaviour of the 
shocks is entirely consistent with the predictions of the classical shock theory. In §|| we 
consider various aspects of the problem of intermediate shocks and discuss ways in which 
to avoid their appearance in MHD simulations with planar symmetry. In particular, 
we present the results of one dimensional simulations using a modified Glimm scheme 
(Glimm 1965) in which these shocks do not appear. 



2. General Theory of Shocks 

In this section we give a brief review of the classical theory of discontinuous solutions 
of hyperbolic conservation laws. For our purposes it is sufficient to consider only the 
dimensional equations of the form 

9u 9f „ . „ . 

^ + ^ = °< ^ 
where u £ lZ n is a vector of conserved variables and f(u) 6 lZ n is a vector of the 
corresponding fluxes. 

As is well known, the system (2T) is called hyperbolic if the Jacobian matrix 

du 

has n real eigenvalues, Xk (k = 1 . . . n) corresponding to n linearly independent right 
eigenvectors, and is called strictly hyperbolic if all the Xk are different. The physical 
significance of the Xk is that they are the speeds of small amplitude waves. 
Waves are classified as linear or nonlinear according to the behaviour of 

Cfc(u) ee r fc (u) • V„A fe (u). 

If Cfc(u) = for all u, then the k-wave is called linear, whereas if the dimension of 
the surface defined by Cfc(u) — is less then n, then it is called nonlinear or genuinely 
nonlinear. 
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The states, u/, u r on either side of a discontinuity travelling with speed s must satisfy 
the shock equations 

s(u { - u r ) = f/ - f r , (2.2) 

The number n s of independent shock equations can be less then n. For example, a contact 
discontinuity in gas dynamics has n s = 3 whereas n = 5. Since A is the Jacobian, we 
clearly have s — > Xk for some k as u; > u r , which means that one can associate each 
discontinuity that allows this limit with one of the waves of the system. A discontinuity is 
called linear if the corresponding characteristic speed does not change across it, otherwise 



it is called nonlinear. There mere fact that a discontinuity satisfies (2.2) it does not 



necessarily imply that it is either stable or that it can arise from continuous initial data. 



For some hyperbolic systems equations (2.2) allow nonlinear shocks that propagate 



with a characteristic speed associated with a nonlinear wave, which means that they can 
be attached to such a wave to form compound waves. Systems with such shock solutions 
are called non-convex. Compound waves may arise from continuous initial data if the 
system allows single simple waves in which Cfe(u) changes sign along the phase curve 
of a simple wave. This condition is therefore often used as an alternative definition of 
non-convexity. Although these definitions are equivalent for a single conservation law, 
they are not necessarily so for systems. 

The evolutionary condition is directly related to the question of existence and unique- 
ness of discontinuous solutions. It is well known that for hyperbolic equations there is a 
general way of deciding this question, which is to use the compatibility conditions that 
must be satisfied along the characteristics (Friedrichs 1955). If a characteristic with wave 
speed Afe enters one side of a discontinuity, then the state on that side must satisfy the 
compatibility relation associated with that characteristic, 

l fe (u)-du = 0, 

where lfc(u) is the left eigenvector of A corresponding to that characteristic. These equa- 
tions are independent provided the are linearly independent i.e. for all hyperbolic 
systems. If the wave speeds on either side of the discontinuity are such that to, com- 
patibility relations have to be satisfied, then there are n s + mi equations relating the 
2n + 1 unknowns associated with the discontinuity, u/, u r and the shock speed, s. A 
discontinuous solution can therefore only exist and be unique if 

mi — 2n — n s + 1. (2-3) 



Obviously, when n s = n, (2.3) reduces to 

rrii = n+l. (2.4) 

It is clear from this that if a characteristic is parallel to the shock curve, then it is 
counted as incoming since the corresponding compatibility relation must be satisfied 
(Gelfand 1963). 

If rrii > 2n — n s + 1 then the system is overdetermined and there is no solution except 
for certain special initial conditions. There will therefore always be arbitrarily small 
perturbations of this data that will destroy such a discontinuity by splitting it into a 
number of waves, just as an arbitrary initial dicontinuity splits in a Riemann problem. 
If mi < In — n s + 1, then the solution exists, but is not unique and one might hope that 
this nonuniqucness can be removed by including dissipative terms. In the following we 
will call condition ( |2.3| ) the strong evolutionary condition and call the condition 

m* < 2n — n„ + 1, 
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which allows nonuniquc solutions a relaxed evolutionary condition. 



An equivalent way of obtaining (2.3) is by a linear structural stability analysis of shock 
solutions (e.g. Landau & Lifshitz 1959, Jeffrey & Taniuti 1964). A discontinuity that 
is exposed to a small amplitude incident wave will only survive if it can respond by 
changing its speed and emitting small amplitude waves. Each such wave is described 
by one parameter and we also have the perturbation in the shock speed, which means 
that there are m + 1 unknowns in this problem, where m is the number of outgoing 
characteristics. Since these are related to the amplitude of the incoming wave by the n s 
shock relations, the discontinuity can only have a unique response if 

m = n s — 1. (2-5) 

It is worth pointing out that, contrary to what is claimed in Myong & Roe, 1997a, this 
analysis does not assume that the discontinuity is weak. This suggests that non-unique 
discontinuous solutions should spontaneously self-destruct by emitting waves even if they 



are not perturbed (Anderson 1963). Although the conditions ( [2.3|) and (2.1) appear 
to be different, the fact that m +m; = 2n means that they are entirely equivalent 
(Gel'fand 1963). Note that, if the system of shock and compatibility equations splits 
into independent subsets, then the discontinuity is only evolutionary if each of these 
subsets has the same number of equations as variables (Jeffrey & Tanuiti 1964). 

Finally, as far as the evolutionary conditions are concerned it does not matter whether, 



or not, the system (2.1) is strictly hyperbolic and convex since these properties are not 



used in the derivation of (2.3,2.5). However, it is only in the case of strictly hyperbolic 



systems that these conditions reduce to the Lax conditions (Lax 1957) 

A fe _i(u ; ) < s < A fe (u ; ) 
A fc (u r ) < s < A fc+ i(u r ) ' 

for a nonlinear discontinuity associated with the fcth characteristic (here we have assumed 
that Ai < A 2 < ... < A n ). 



3. Evolutionary conditions and dissipative shock structure 

In order to assess recent claims that nonevolutionary shocks become admissible if dis- 
sipative terms are included, we need to look at the general relationship between the 
evolutionary conditions and the uniqueness and existence of steady dissipative shock 
structures. Godunov (1961) has shown that it is much easier to explore this question if 
the equations can be transformed to a symmetric form. Although this is not possible 
for arbitrary hyperbolic systems of conservation laws, it can certainly be done for gas- 
dynamics, MHD, and the shallow water equations and probably for any system that can 
arise in nature. 

3.1. Symmetric Form of the Ideal Equations 

We start by summarizing some of the results described by Friedrichs (1954), Friedrichs 
& Lax (1971) and Boillat (1974, 1982). As before, it is only necessary to consider the 
one dimensional case. 

Consider a dissipation-free system of conservation laws described by the equations 
( |2.ip . Suppose now that there exists a quantity, h(u), which is also conserved as long as 
the solution to this system is continuous. For example, h(u) is the entropy in gasdynamics 
or MHD, whereas it is the total energy for the shallow water equations. If such a quantity 
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exists, then there must exist a flux function, g(u), such that 



( |2.l| ) and ( |3 . l|) can only be consistent if 

dh_djj_ _ dg 

dui duj duj ' 



(summation convention assumed), since then 

dh dg_ _ dh_ ( chn df^ _ 
dt dx dui \ dt dx 



In terms of the variables u', (|2.l|) becomes a symmetric system 

dt dx 

where the symmetric matrices P and Q are given by 

dui d 2 ti d 2 h 



du'j du'idu'j duiduj 



Qi 



dh d 2 g' 



(3.2) 



(3.3) 
(3.4) 



for any C 1 solution satisfying (2.1). 

If we now use h to define the Legendre transformation 

' - dh 

dh 1 

U * = ^> 

h' = h + u'iUi, (3-5) 
then ( |3.2j ) allows us to write the fluxes as 

~ d< 

where 

g' = g + u'ifi- 



(3.7) 



du'j du\du'j 

Note that h is usually a strictly concave function, in which case (3/7) ensures that P 
is positive definite and the transformation is non-singular. In ordinary gasdynamics or 
MHD, h is the entropy per unit volume and is therefore guaranteed to be concave by the 
second law of thermodynamics. For the shallow water equations h = — e, where e is the 
sum of the kinetic and potential energy and dissipation ensures that this is also concave. 

3.2. Dissipative Equations 
If we now assume that the dissipative fluxes are proportional to the spatial gradients of 



the dependent variables, then the dissipative version of (3.6) is 
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where D is a matrix of dissipation coefficients. Multiplying this on the left by u" (the 
superfix t denotes the transpose) and using (3.1-3.3) gives the evolution equation for h 

dh dg , t d j^du' 

dt dx dx dx ' 

Integrating this over an arbitrary fixed interval [a, b] and integrating the dissipative 
term by parts gives 



— / hdx 

dt 



u'*D 



du' 

dx 



<9u'* <9u' , 
dx - 

ox ox 



Since the term on the RHS of this equation represents a source term for h and the second 
law of thermodynamic requires that this be positive if h is the entropy per unit volume, 
the matrix D must be positive definite for gasdynamics and MHD. The dissipative shallow 
water equations must also satisfy this condition if we set h — — e, where e is the total 
energy. 

One can also show that all lin ear waves decay if D is positive definite and h is a strictly 
concave. The linear version of (3^) is simply 

p<9u' ^ q^ u ' d^ 2u ' 

dt dx dx 2 

where P, Q, and D are now constant matrices. Multiplying this by u" and integrating 
over [a, b] gives 



4 [ u'M'u 

dt 



du' 

u"Qu' 2u"D^- 

dx 



b du n du' 

dx dx 



after integrating the dissipative term by parts. Since P is positive definite if h is strictly 
concave, the term on the RHS ensures that all linear waves decay if D is positive definite, 

3.3. Steady Shock Structures 
Now consider a solution of the steady version of fl3.8p 



dx 



d d , 

dx dx 



with the boundary conditions 



u 



u' x — > 



-00 

\-oo. 



(3.9) 



(3.10) 



If this represents a shock structure, then u[ and u^, must satisfy the shock relations in 
the shock frame 



f K) = f «). 



Integrating (3.9) and applying the boundary conditions (3.10) gives 

du' 



D- 



dx 



f(u')-f(uD=f(u')-f«) 



(3.11) 



(3.12) 



A steady shock structure therefore corresponds to a solution of (3.12) that connects the 
equilibrium points uj and u^.. We now show that there is no guarantee that this solution 
is unique and structurally stable unless the corresponding discontinuous solution of the 
ideal system are satisfies the evolutionary conditions (2.3). 

Let L u be the unstable manifold of the point uj and R s the stable manifold of the point 
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u' r . Then the trajectories in L u and R s are described by dim(L u ) — 1 and dim(R s ) — 1 
parameters respectively. Since any trajectory which lies in both has to satisfy n — 1 
matching conditions, this means that, in general, there will only be a unique trajectory 
connecting uj and u r if dim(L u ) + dim(R s ) = n + 1. If dim(L u ) + dim(R s ) > n + 1, 
then the trajectory may not be unique, whereas if dim(L u ) + dim(R s ) < n + 1, then 
any trajectory that does exist can be destroyed by perturbations of uj, u^. i.e. it is not 
structurally stable. 

The following theorem relates dim(L u ) and dim(R s ) to the number of characteristics 
entering the shock: 



Theorem 3.1. //uj, is an equilibrium point of the dissipative shock equations (3.1i) 
at which none of the characteristic speeds vanish, then the equilibrium point is hyperbolic 
and the dimension of its stable ( unstable ) manifold is given by the number of positive 
(negative) characteristic speeds in the state u' e . 

Proof. Suppose that u' e = uj (the proof for u^. is identical). Then linearizing (3.12) 
in the neighbourhood of u[ gives 



D ; — = Q ; v, 

ax 



where v = u' — uj, Q; = Q(uJ) and D; = D(uJ). If this equilibrium point is hyper- 
bolic, then the dimension of its stable (unstable) manifold are given by the numbers of 
eigenvalues, satisfying 

|Q;- M D;| = 0. (3.13) 

with positive (negative) real parts. 



On the other hand, the characteristic speeds for the system (3.6), A^, in the state u[ 
are given by 

|Q,-AP ; |=0. (3.14) 

A standard result (e.g. Gantmacher 1959) tells us that, since P;, Q; are symmetric and 
Pi is positive definite, Q; has the same number of positive, negative and zero eigenvalues 
as the set A^. If, like Godunov (1961), we assume that D; is symmetric as well as positive 



definite, then the theorem would follow immediately from (|3.13 ) and ( 3.14 ). However 



the following lemma shows that this is an unnecessary restriction. 

Lemma 1. Let Q be a non-singular symmetric matrix, D a positive definite matrix 
and fik the solutions of 

|Q-/xD| = 0. 

Then the number of with positive ( negative ) real part is the same as the number of 
positive (negative) eigenvalues o/Q. 

Proof. Define 

D e =D s +eD Q , 

where eg [0,1] and 

D s = i(D + D*), D Q = i(D-D*). 

It easy to see that D e is also positive definite. 
Now consider the eigenvalue problem 

|Q-/i(e)D e |=0. 
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The conclusion of the lemma is certainly true for e = 0, since then D e is symmetric. If 
we can show that the /^&(e) are continuous functions of e and that 9?{/^fe(e)} ^ Vfc for 
e e [0, 1], then it will also be true for e = 1. 

The Hk (e) are the roots of a polynomial of degree n whose coefficients are polynomials 
in e. A root can therefore only change discontinuously by going to infinity, which can 
only occur if the coefficient, |D e |, of the highest power of u vanishes. However, this 
cannot happen since D t is positive definite for e € [0, 1]. The /Ufc(e) must therefore be 
continuous functions of e for e € [0,1]. 

In order to prove that the fik cannot cross the imaginary axis, suppose that for some 
k, /Ufe(e) = if], where rj is real. If a + ih is the corresponding eigenvector, we have 

Qa + r/D £ b = 0, 
Qb-?7D e a = 0. 

Multiplying the first of these by b* . the second by a* and substracting gives 

^(b'D.b + a*D e a) = 0. 

Since D e is positive definite this requires rj = and hence ^ — 0, which cannot be true 
if the eigenvalues of Q are non-zero. This completes the proof of the lemma. □ 

( |3.13 ), flS.14 ) and lemma (Q) show that the theorem is true even if D is not symmetric. 

□ 

This is a somewhat more direct proof of a result which has also been obtained by 
Kulikovsky & Lyubimov (1965). In their analysis of viscous shock structures, Myong & 
Roe (1997a) assumed that Theorem 3.1 holds for MHD, but did not give a proof. 



This analysis tells us that if the shock relations (3.11) have a solution such that none 



of the characteristic speeds given by (3.14) vanish in both the left and the right state 



and rrii is the number of characteristics entering the shock, then 

(a) for m-i = n+1 the shock can have a unique structurally stable dissipative structure; 

(&) for m% > n + 1 the dissipative structure is not guaranteed to be unique. 

(c) for mi < n + 1 there might be a unique dissipative structure but it cannot be 
structurally stable. 

These conditions are not only compatible with the evolutionary conditions, they are 
complementary to them. Shocks for which rm > n + 1 have a dissipative shock structure 
and could therefore be regarded as admissible on these grounds. However, the left and 
right states of such shocks must be carefully tuned since they cannot adjust themselves 
to an arbitrary small perturbation of their left and right states. Shocks that satisfy the 
relaxed evolutionary condition, < n + 1, are apparently permitted by the ideal equa- 
tions, but cannot establish a dissipative structure and must spontaneously self-destruct. 
It is therefore clear that the only physically admissible shocks are those those that satisfy 



the strong evolutionary conditions (2.3) or (2.4) 



Theorem 3.1 gives us no information in those cases for which the shock speed coincides 
with at least one of the characteristic speeds. The corresponding critical point is then 
no longer hyperbolic and its type depends on the details of the particular system. 



4. Application to Magnet ohydro dynamics 

As we shall see, the mathematical properties of the full system of MHD and the reduced 
planar system of MHD are somewhat different and this has to be clearly understood when 
the evolutionary conditions are applied. We therefore discuss these systems separately. 
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4.1. Full System of MHD 

It is well known that the one dimensional equations of MHD can be written in the form 
( |2.lD (e.g. Brio & Wu 1988). The conserved quantities u and the corresponding fluxes f 
are 



u 



p 




pv x 


PV X 




pvl+ Vg + B 2 /2-Bl 


P Vy 




pV X Vy ~ B X By 


pv z 


f = 


pv x v z - B X B Z 


e 




{e + p g + B 2 /2}v x -B x (v.B) 


By 




V X By - VyB X 


. B * _ 




v x B z - v z B x 



Here p g is the gas pressure, 



1 9 

2 B 



1 2 

- 2 pv 



is the total energy per unit volume and i is the enthalpy per unit volume. Here we use 
units such that the velocity of light and the factor Air do not appear. 

As we have already discussed, ideal MHD has a supplementary conservation law repre- 
senting the conservation of thermodynamic entropy. The second law of thermodynamics 
guarantees that the function h = pS, where S is the entropy per unit mass, is strictly 
concave (e.g. TerHaar & Wergeland 1966) and hence that the matrix P defined by (3.7) 
is positive definite. The system of MHD equations can therefore be written in the sym- 



metric form (3.6) and is hyperbolic. Although this has been demonstrated for relativistic 
MHD by Ruggeri & Strumia (1981), we have been unable to find an account of the cor- 
responding analysis for classical MHD in the literature. However, since the derivations 
are similar to those for the relativistic case, we shall simply give the symmetric variables. 
They are 



1 
T 



1 



T ' 
B. 



T ' 



T ' 



1 i B y ! B z 

There is no need to verify that the matrix, D, of dissipation coefficients is positive 
definite, since this must be true for any system that obeys the second law of thermo- 
dynamics. Indeed, this condition is used to derive the dissipative equations in the first 
place (e.g. Landau & Lifshitz 1960). The exact form of symmetrized equations is also of 
no importance for our purposes. Their existence, does, however, allow us to apply the 
conclusions of the general theory described in Sections 2 and 3 to dissipative MHD. 

4.1.1. Characteristic Wave Speeds 

Since there are seven variables in this system, there are seven waves whose speeds are 



Fast Waves 
Alfven Waves 
Slow Waves 
Entropy Wave 



v x T cf, 

V x T C a , 

v x T c s , 
v~. 



where the alfven speed, c a , and the slow and fast speeds, c s , c/ are given by 

c a = \B x \/y/p, 
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B*\ 2 Aa 2 B 2 S 
P 



2 1 



2 




1/2- 



where a is the adiabatic sound speed. Note that < c s < c a < Cf. If B x = then 
c s = c a = 0, whereas if the transverse component of the magnetic field, Bt, vanishes 
then Cf = c a if c a > a, c s = c a if c a < a and c s ~ Cf = c a if c a = a. The MHD 
equations are therefore not strictly hyperbolic. Brio & Wu (1988) also argued that they 
are non-convex, but we shall postpone discussion of this until later. 

4.1.2. Shock Types 

The MHD shock equations allow two linear solutions and several distinct types of 
nonlinear solutions which satisfy the entropy principle that the entropy of a fluid element 
always increases. A convenient way of classifying these is to use the jump in the transverse 
component of magnetic held, B t . From the shock equations one finds (Jeffrey & Taniuti 
1964) 

[B t (c 2 a -v 2 x )] l = [B t (cl-v 2 x )] r , (4.1) 

where v x is the velocity in the shock frame. Note that if c 2 — v 2 does not vanish, then 
Bt on one side of the discontinuity must be either parallel or anti-parallel to that on the 
other. 

The nonlinear solutions are 

(a) Slow/Fast shocks, which have non-zero Bt in the same direction on both sides. 



(4.1) then implies that there is no change in sign of (c 2 — v 2 ). The magnitude of magnetic 
field is larger on the downstream side for fast shocks and smaller downstream for slow 
shocks. 

(b) Intermediate shocks, which also have non-zero B t but in opposite directions on 
either side of the shock (Anderson 1963, Cabannes 1970). ( |4.1| ) then implies that (c 2 — v 2 ) 
changes sign. 



(c) Switch-on shocks, which have vanishing B t upstream. (4.1) then implies that 
I = c 2 on the downstream side. 



(d) Switch-off shocks, which have vanishing B t downstream. (4.1) then implies that 
v 2 = c 2 on the upstream side. 
The linear discontinuities are 

(a) Alfven discontinuities, which have v 2 = c 2 on both sides, case (4.1) then allows 



an arbitrary change in the direction of B t . However, the magnitude of B t remains 
unchanged, which is why these are sometimes called rotational discontinuities. 

(&) Contact discontinuities, which have the same value of v x on both sides, but v 2 ^ c 2 . 



(4.1) then requires that B t be continuous unless B x — and the other shock conditions 
require all other variables, except for the density, to be continuous. 

We shall also find occasion to use the following classification of nonlinear MHD shocks, 
which is due to Germain (1960). The states in the shock frame are divided into four types 

1) \v x \ > Cf, 2) c f > \v x \ > c a , 

3) c a > \v x \ > c s , 4) c s > \v x \, 

and a shock is defined to be of type m — > n if the upstream and downstream states are 
of types m and n respectively From the MHD shock equations, one finds that pressure 
and specific volume, r (r = 1/p), on each side of a nonlinear shock satisfy the following 
equations 
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Figure 1. The shock invariant, H, as a function of specific volume, r, for three different cases. 

2 2t q (t - r a y 

where G is the mass flux, F x , F y , H are shock invariants and r a = B^/G 2 . The analysis 
in Anderson (1963) can be used to show that the function H(t) is as shown in Figure 
1. r — Tj has the same sign as v 2 — cf, where i = s,a, f. One can see that there are 
six different types of compressive shocks: fast shocks (1 — > 2), slow shocks (3 — > 4), 
and four intermediate shocks: 1 — > 3, 1 — > 4, 2 — * 3, and 2 — * 4. Depending on 
the relative position of the maxima of H, there are also limit shocks which propagate 
with the fast speed relative to the upstream state and/or the slow speed relative to the 
downstream state (see figures lb,c). We shall denote such such shocks by / — * n and 
n — > s respectively. These shocks turn out not to be evolutionary, but if they were, then 
MHD would be a non-convex system. 

4.1.3. Evolutionary conditions 

When we apply the evolutionary conditions to MHD discontinuities we have to take 
into account the fact that the system of shock and compatibility equations split into two 
independent subsets for all types of discontinuities, except the alfven discontinuity,. If 
we choose a reference frame such that on one side of a discontinuity B z = and v z = 
then the system of shock equations contains two equations involving B z and v z . These 
are 

B z i = B Zr , 

V z i = V Zr . 

The compatibility relations along the alfven characteristics only involve B z and v z and 
they also the only ones that do so. An evolutionary discontinuity that is not an alfven 



discontinuity must therefore not only satisfy the general condition (2.3), but also have 
exactly two incoming, and hence two outgoing alfven characteristics. These conditions 
also follow from the linear stability analysis (Syrovatskii 1959, Jeffrey & Tanuiti 1964). 

In the rest of this section we simply state the well known results on the evolutionary 
properties of MHD discontinuities. We do, however, pay particular attention to those 
cases in which there are characteristics travelling with the same speed as the discontinuity. 
As we have pointed out in Section 2, such characteristics must be counted as incoming. 

There is no dispute about the fact that fast and slow shocks are evolutionary because 
they have eight incoming characteristics, two of which are alfven waves. Furthermore, 
since their speed can never be equal to a characteristic speed, Theorem 3.1 tells us that 
they also have a unique structurally stable dissipative structure. 

All intermediate shocks are super-alfvenic with respect to the upstream state and sub- 
alfvenic with respect to the downstream state, which means that they have too many 
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(> 2) incoming alfven characteristics. They are therefore nonevolutionary and can be 
destroyed by interactions alfven waves. 

The same argument applies to switch-on and switch-off shocks which also have too 
many (9) incoming characteristics, 3 of which are alfven characteristics. However, these 
solutions are clearly limits of fast and slow shocks and therefore have evolutionary solu- 
tions in their immediate neighbourhood, which is why Jeffrey & Taniuti (1964) call them 
weakly evolutionary. That they are not strictly evolutionary can also be understood from 
the following example. Consider a switch-on shock overtaking weak switch-off fast rar- 
efaction travelling in the same direction. Once these have merged, the shock is no longer 
propagating into a state with zero transverse magnetic field. Since the shock is superfast, 
it has no way of modifying its upstream state and therefore cannot remain a switch-on 
shock. Instead, such an interaction leads to the appearance of a neighbouring fast shock 
solution, together with some other waves, at least one of which must, in general, be an 
alfven wave. 

If we count the two entropy characteristics as incoming on the grounds that they have 
the same speed as the discontinuity, then contact discontinuities have eight incoming 
characteristics, two of which are alfven characteristics. They are therefore evolutionary. 

Alfven discontinuities also have eight incoming characteristics if we include the two 
alfven characteristics that have the same speed as the discontinuity. The total number 
of incoming alfven characteristics is three, but this is allowed since the fact that the 
shock equations for these discontinuities couple the y and z components of velocity and 
magnetic field means that this is the one case for which the shock equations do not 
decompose into two sets. 

Theorem 3.1 cannot be applied to contact and alfven discontinuities since they prop- 
agate with a characteristic speed. However, they would in any case not possess a steady 
dissipative structure simply because they are linear and therefore have no nonlinear 
steepening to balance the spreading due to dissipation. For this reason, Wu (1988b) 
considers them to be inadmissible, but since their width grows like t 1 / 2 , whereas the 
separation between the waves in a Riemann problem grows like t, they must be regarded 
as admissible components of the solution for large times. 



4.2. Reduced system of planar MHD 

In this section we discuss the system of equations which describes MHD in a world in 
which the plane defined by the velocity and the magnetic field is invariant. There are 
several reasons for doing this. Firstly, it has some interesting properties. Secondly, we 
want to show that the general classical theory of shocks is as valid for this system as it 
is for the full system. Finally, the numerical simulations that gave rise to the current 
conroversy surrounding intermediate shocks reflect the properties of this system. 

When the z components of the magnetic field and velocity vanish, the equations reduce 
to a system of 5 variables with the following vectors of conserved quantities and fluxes 



p 




PVx 




P Vy 


f = 


e 




By 





Bl 



{e+Pg 



pv 2 x + Pg + B 2 /2- 

pV X Vy - B X By 

+ B 2 /2}v x -B x (v.B) 

V X By VyB X 



This is still a hyperbolic system but it is fundamentally different from the full system 
of MHD because it does not have alfven waves. However, the other characteristic fields 
are still present with the same eigenvalues and with eigenvectors that are the same 
apart from the reduced number of components. Moreover, it has the same solutions of 
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the shock equations including the alfven discontinuity, except that these are now only 
allowed to change the direction of the transverse magnetic field by ir. This follows from 
the remarkable property of the full system of MHD that there exists an inertial frame 
in which the variation of the transverse components of the magnetic field and velocity 
induced by all characteristic waves and shocks, except for alfven waves, are confined to 
single plane. Note that the alfven discontinuity still propagates with the alfven speed, 
but this is no longer on of the characteristic speeds. The Ricmann problem for this 
system has been analysed in considerable detail by Myong & Roe (1997b) who came to 
the conclusion that the classical evolutionary conditions are inadequate for this system. 
However, we intend to show that this claim is based on a failure to recognise the essential 
difference between the reduced system and full MHD. 

4.2.1. Evolutionary conditions 

Since the number of equation is reduced by two and it is the alfven waves that are 
lost, we can conclude that all evolutionary discontinuities that have two incoming alfven 
characteristics in the full system remain evolutionary in the planar system. This implies 
that fast, slow and contact discontinuities are evolutionary. 

On the other hand, discontinuities that are evolutionary in the full system, but which 
do not have exactly two incoming alfven characteristics must be non-evolutionary in the 
planar system. There is only one such discontinuity, the alfven discontinuity, which now 
only has 5 incoming characteristics and should therefore spontaneously self-destruct even 
if it is not perturbed. 

Another interesting feature is that some of the shocks that are non-evolutionary in the 
full system become evolutionary in the reduced system. 1 — > 3 shocks now satisfy the 
strong evolutionary condition, in fact they have the same incoming and outgoing char- 
acteristics as fast and switch-on shocks. As far as the characteristic count is concerned 
these three shocks arc therefore indistinguishable so that one can use a single name, 
plane fast shock, say, for all of them. Similarly, 2^4 shocks, switch-off shocks and slow 
shocks become slightly different versions of evolutionary plane slow shocks. 

However, 1 — ► 4 shocks remain non-evolutionary even in the plane system since they 
have 7 incoming characteristics. Such shocks, which have too many incoming charac- 
teristics, are often called over compressive in the literature. As we have shown, although 
they do have a steady dissipative structure, it is not unique and it does not help them 
to survive interactions with external perturbations. 

2^3 shocks have only 5 incoming characteristics and are therefore non-evolutionary. 
Such shocks, which have too few incoming characteristics, are often called undercom- 
pressive. Since they do not have a structurally stable steady dissipative structure they 
should disintegrate spontaneously even without any external perturbation. 

Now consider shocks that propagate at one of the characteristic speeds in either the up- 
stream or downstream state. 1 — > s, f — > 4, and / — ► s shocks are non-evolutionary since 
they have 7 incoming characteristics. On the other hand, 2 — * s and / — ► 3 shocks have 
6 incoming characteristics and are therefore evolutionary. The planar system of MHD 
is therefore genuinely non-convex and admits the two evolutionary compound waves: a 
slow compound wave consisting of a 2 — ► s shock with an attached slow rarefaction and 
a fast compound wave consisting of a fast rarefaction with an attached / — > 3 shock. 

Finally, we list the evolutionary shocks and compound waves of the planar system 
along with the notation used in Myong & Roe (1997b): 

Slow planar shock (SI); 

Fast planar shock (S2); 

Slow compound wave (CI); 
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Fast compound wave (C2); 

Contact discontinuity (not considered). 
Myong & Roe (1997b) found that some Riemann problems only have a solution if non- 
evolutionary shocks are permitted. However, as we discuss in these Riemann problems 
are confined to regions of parameter space with zero volume, which is exactly what is 
meant by the statement that non-evolutionary shocks are structurally unstable. 

In the next section we show that the results of numerical calculations are entirely 
consistent with these conclusions. 



5. Numerical Calculations 

The numerical calculations were carried out using the scheme described in Falle, Komis- 
sarov & Joarder (1998). This is an upwind shock capturing scheme which is capable of 
dealing with shocks of arbitrary strength even without the inclusion of any dissipation 
other than that introduced by the truncation errors. Careful test simulations have shown 
that this scheme provides accurate solutions for all types of MHD waves in all regimes. 
One can argue that if a numerical scheme works well then its numerical dissipation must 
have the same qualitative properties as the physical dissipation. However, in order to 
remove any doubts, we modified our scheme so that it can now handle dissipative MHD 
and all the calculations described here have a fully resolved dissipative shock structures 
(about 15 mesh poin ts w ide). For this we used a simple scalar form for the dissipation 
for which equations ( |2.l|) become 

du df d 

dt dx dx^ 



where the diffusive fluxes are 
/ 





4^ dv x 
3 dx 



\ 



/' 



A[iv x dv x 
3 dx 



d_Uy 

dx 



dVy 

dx 

dv z 
U dx 

dv z 
dx 

dB v 



By -dx~ 



B 7 



dB z 
dx 



dx 
dB z 

V *»ftT J 

where fj, is the dynamic viscosity, k the thermal conductivity and v m the resistivity. 

As expected, the outcomes of all the simulations presented here did not not depend 
on the size of dissipation and were the same even when only numerical and/or artificial 
dissipation was present. The only effect of changing the dissipation was to alter the form 
and width of the shock structures. 

First of all, we need to establish whether the behaviour of numerical MHD shocks 
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2^3 Intermediate Shock (figure 2a) 
Left state: p = 1, p g = 1, v = (-0.95, 0, 0), B = (1, 0.5, 0) 
Right state: p = 0.837, p g = 0.705, v = (-1.135, 1.266, 0), B = (1, -0.7, 0) 

Alfven Shock (figure 2b) 
Left state: p = 1, p g = 1, v = (-1, 1, 0), B = (1, 1, 0) 
Right state: p = 1, p 9 = 1, v = (-1, 3, 0), B = (1, -1, 0) 

1 — » 3 Intermediate Shock (figures 3a, 5 and 7a) 
Left state: p = 1, p g = 1, v = (-0.925, 0, 0), B = (1, 0.5, 0) 
Right state: p = 0.498, p g = 0.258, v = (-1.857, 0.648, 0), B = (1, -0.1, 0) 

2^4 Intermediate Shock (figures 3b and 7b ) 
Left state: p = 1, p g = 1, v = (-0.4, 0, 0), B = (0.5, 0.5, 0) 
Right state: p = 0.561, p g = 0.155, v = (-0.714, 2.252, 0), B = (0.5, -1.3, 0) 

1 — » 4 Intermediate Shock (figure 4) 
Left state: p = 1, p g = 1.2, v = (-0.842, 0.0, 0.0), B = (1.0, 0.4, 0) 
Right state: p = 0.390, p g = 0.161, v = (-2.16, 0.644, 0), B = (1.0, -0.142, 0) 

Brio & Wu Problem (figure 8) 
Left state: p = 1, p g = 1, v = (0, 0, 0), B = (0.75, 1, 0) 
Right state: p = 0.125, p g = 0.1, v = (0, 0, 0), B = (0.75, -1, 0) 



Table 1. Riemann problems for the numerical calculations. 



Problem Domain n p/p n/p v, 



Figure 2a 


h4,l] 


Figure 2b 


[-2,1] 


Figure 3a,b 


[-4, 1] 


Figure 4 


[-4, 1] 


Figure 5 


[-1,1] 


Figure 6 


[-2,1] 


Figure 7a 


[-8,2] 


Figure 7a 


[-14, 1] 


Figure 8a,b 


[2.5,4.5] 



250 0.02 0.01 0.01 

150 0.02 0.01 0.01 

250 0.02 0.01 0.01 

250 0.02 0.01 0.01 

200 0.01 0.005 0.005 

300 0.01 0.005 0.005 

500 0.02 0.01 0.01 

750 0.02 0.01 0.01 

200 0.0 0.0 0.0 



Table 2. Other parameters for the numerical calculations, n is the number of mesh points, p 
is the kinematic viscosity, k is the thermal conductivity, v m is the resistivity. 



agrees with the predictions of the evolutionary theory. In order to do this, we adopt the 
following procedure. First we test whether a shock has a steady dissipative structure 
by setting up the relevant Riemann problem and running the calculation until a well 
resolved steady dissipative shock structure is established, as expected for evolutionary 
and overdetermined shocks, or a completely different solution emerges, as expected for 
underdetermined shocks. If a steady structure exists, then we test to see whether it can 
survive small perturbations. This can be accomplished by considering a slightly different 
Riemann problem, as in Barmin ct al. (1996) or, like Wu (1988a), allowing a small 
amplitude wave to interact with the shock. 
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Figure 2. Planar simulations of shocks that should not have a steady dissipative structure in 
planar MHD: (a) 2 — > 3 shock, (b) AlfVen shock. In both cases the outcome is a slow compound 
wave (SCW). The dashed lines show the corresponding initial solutions. The continuous lines 
show the final solutions. 



5.1. Planar MHD 

We start by discussing the results of the planar simulations. They show that if the initial 
discontinuity corresponds to a slow planar shock then a smooth steady shock structure 
connecting the initial left and right states finally develops and it does not matter whether 
the shock is 3 — > 4 or 2 — > 4. The same thing happens for the fast planar shock and the 
overdetermined (overcompressive) 1 — ► 4 shock. In contrast, figure 2 shows that 2^3 
shocks and alfven shocks always turn into a slow compound wave. All this is exactly 
as perdicted by the theory described in ^ and Our simulations cannot be used to 
determine whether limit shocks, such 1 — > s and / — > 3), have a steady dissipative shock 
structure, simply because it is impossible to set up a shock whose speed is exactly equal 
to a characteristic speed. However, if we compute a Riemann problem that corresponds 
to a compound wave of any of the types discussed above, the wave that is expected, 
or strictly speaking a solution close to such a wave, always emerges. This is hardly 
surprising because all of them have neighbouring solutions containing shocks with a 
steady dissipative structure. 

As is shown in figure 3, evolutionary shocks always survive interactions with small 
amplitude waves and persist if the Riemann problem is perturbed. Figure 4 shows how 
a small variation of the initial data forces an overdetermined 1 — > 4 shock to split into 
two evolutionary shocks. Depending on the form of the perturbation, the shock either 
splits into a 1 — > 2 shock followed by a 2 — ► 4 shock or a 1 — > 3 shock followed by a 3 — > 4 
shock. This is to be expected because, as one can see from figure 1, a 1 — > 4 shock is 
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Figure 3. Planar simulations of the interaction between evolutionary shocks and small ampli- 
tude fast rarefactions (SB t = 10%). (a) fast (1 — ► 3) shock, (b) slow (2 — > 4) shock. In both 
cases the outcome is a shock of the same type, together with some other waves. Here FR denotes 
a fast rarefaction and V x is the x-component of velocity as measured in the shock frame. The 
dashed lines show the initial solutions. The continuous and dotted lines show the final solutions. 



exactly equivalent to one or other of these shock pairs propagating with the same speed. 
In fact, this result is in complete agreement with the analysis of the Riemann problem 
for planar MHD in Myong and Roe (1997b). 1 — > 4 shocks, O-shocks in their notation, 
are only required on the boundary between the two domains of parameter space in which 
their solution involves a combination of fast and slow planar shocks (S2 and SI). 

The results for compound waves involving non-evolutionary shocks are similar. Fig- 
ure 1 shows that the non-evolutionary 1 — > s limit shock can be understood as a double- 
layer shock composed of two evolutionary shocks, a 1 — > 2 and a2^s. Indeed, if the 
Riemann problem corresponding to a compound wave containing such a shock is per- 
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Figure 4. Planar simulations of 1 — ► 4 shock subjected to a small variation of pressure (±10%) 
in the left state. This shock is non-evolutionary even in planar MHD and splits as the result of 
the perturbation into two evolutionary shocks plus other small amplitude waves The outcome is 
(a) 1 -> 3 and 3 -> 4 shocks if Sp = -10% and (b) 1 -» 2 and 2 -» 4 shocks if Sp = +10%. The 
dashed lines show the initial solutions. The continuous and dotted lines show the final solutions. 
V x is the x-component of velocity as measured in the frame of the emerged intemediate shock. 



turbed, then in some cases the outcome is a 1 — * 2 shock and a slow compound wave and 
in other cases it is a 1 — > 3 shock and a detached slow rarefaction. 

All this can be summed up by saying that in planar MHD the behavior of shocks in 
our numerical simulations is entirely consistent with the classical evolutionary theory of 
shocks and the theory of dissipative shock structures as described in §[| and §H. 
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5.2. Full MHD 



Since both fast (1 — > 2) and slow (3 — > 4) shocks satisfy the strong evolutionary condition 
in full MHD they are expected to have unique dissipative structure and be stable with 
respect to small perturbations of any kind. This is precisely what we find from our 
simulations. 

1 — ► 3 and 2 — * 4 shocks are overdetermined in full MHD and it is therefore possible 
that they might have a nonunique steady dissipative structure, indeed it turns out that 
they do. These shocks, as well as 1 — ► 4 shocks, can now have a nonvanishing z-component 
of magnetic field inside the shock layer even if B z = outside. For given the dissipative 
coefficients their stucture can be parameterised by the value of the following integral 



We can gradually increase or decrease the value of I z by sending from the downstream 
side of the shock an alfven wave that first rotates the magnetic field by a small angle and 
then restores the original state. This wave is absorbed by the shock which develops a 
new steady structure (see Figure 5a). However, like Kennel et al. (1990) we found that 
there is a maximum value of \I Z \ that the shock can manage. If this limit is exceeded, 
then the shock disintegrates (see Figure b). This does not occur in the case of fast and 
slow shocks because the alfven waves do not get trapped inside the shocks, but instead 
pass straight through. 

2^3 shocks have the right number of incoming characteristics and may therefore 
have a unique dissipative structure in full MHD. Since such a structure does not exist 
in planar MHD, we can only expect to find them in our simulations by allowing a non- 
zero B z . In order to do this, we modified the initial data by inserting a layer in which 
the transverse field rotates smoothly from that in the original left state to that in the 
original right state. We found that the solution never relaxed to a smooth steady 2^3 
transition and were about to conclude that no steady structure exists until we realised 
that the solution shown in figure 56 actually contains a 2 — > 3 shock, which was produced 
by the disintegration of the 1 — > 3 shock. We therefore we studied the reaction of a 1 — > 3 
shock to an increase in I z . After absorbng another alfven wave the shock splits and one 
of the emerging waves is again a 2 — > 3 shock but of smaller amplitude (figure 6). This 
behaviour is consistent with the existence of a unique dissipative structure for 2 — > 3 
shocks. In fact, what happens is that, as I z increases, the shock tends to an alfven shock 
that rotates the transverse field by 7r 

Finally, we have also verified that all intermediate shocks and compound waves disin- 
tegrate when exposed to perturbations that render the left and right states non-coplanar. 
For example, figure 7 shows how 1 — > 3 and 2^4 shocks split into evolutionary waves 
after interaction with a small amplitude alfven wave. After the alfven wave has been 
absorbed the transverse fields on either side of the shock are no longer parallel or an- 
tiparallel as required by the shock equations. The shock can only become coplanar by 
emitting alfven waves, which, for an intermediate shock, can only be done in the down- 
stream direction. However, since there is no downstream travelling alfven wave that can 
restore the original post-shock state, the shock must split. This argument is not new, in 
fact it was used by Kantrowitz & Petschek (1966) to prove that intermediate shocks are 
unphysical. The wave designated as AW in figure 7 can be called a dissipative alfven wave 
but it could also be described as an evolving 2^3 shock with a gradually increasing 
value of I z . 

We therefore conclude that for full MHD the behaviour of shocks in our numerical 
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Figure 5. Dissipative structure of a 1 — > 3 shock in full MHD for different values of I z . This 
shock has a nonunique steady dissipative structure that depends upon J z . For relatively small 
values of I z this structure is steady (left panel, I z = —0.085) but for larger values it splits into 
1 — » 2 and 1 — * 2 shocks (right panel, I z — —0.20. 



simulations is also entirely consistent with the classical evolutionary theory of shocks ana 
the theory of dissipative shock structures as described in §[| and §fx 



6. Discussion 

The results described in the previous sections have clarified many aspects of the shock 
theory in general and MHD shocks in particular and provide a basis upon which we can 
discuss other important, related, issues. 

6.1. Riemann problems and evolutionary conditions 

One of the arguments in favour of non-evolutionary shocks used in current literature is 
that some Riemann problems do not have a solution unless non-evolutionary shocks are 
admitted (e.g. Glimm 1988, Myong & Roe 1997a, b). This is presumably based on the 
belief that any Riemann problem must have a physically admissible solution. Although 
this is certainly true for gas dynamics, there is surely no reason why this has to hold 
for any system. It all comes down to the notions of bifurcations and structural stability. 
One has to ask the following question: is it, or is it not, possible to carry out the relevant 
experiment in a laboratory? If the qualitative result of the experiment does not change 
when the initial conditions are slightly changed, then the problem is structurally stable 
and the experiment is possible, at least in principle. However, if this is not true, then the 
problem is structurally unstable and no appropriate experiment is possible. It therefore 
follows that the set of structurally unstable Riemann problem are confined to regions of 
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Figure 6. Dissipative structure of a 2 — + 3 shock in full MHD. This shock has a unique steady 
dissipative structure and therefore reacts to a change in I z by emitting some waves and turning 
into a different 2^3 shock. The continuous lines show the solution for 7 Z = —0.20 and the 
dashed lines for I z — —0.52. 



parameter space whose total volume is zero. Now suppose there is an MHD Ricmann 
problem that has no other solutions than those containing non-evolutionary shocks. Since 
there are arbitrary small perturbations of the parameters that cause these shocks to split 
into evolutionary shocks, this Riemann problem must be structurally unstable. In full 
MHD the only known case for which a non-evolutionary shock, a 1 — * 4 shock, is required 
is a piston problem in which the piston velocity is parallel to the magnetic field (Jeffrey 
& Taniuti 1964). If this condition is not exactly satisfied then the non-evolutionary shock 
does not arise. Close inspection of the solution of the Riemann problem for planar MHD 
presented by Myong & Roe (1997b) shows that non-evolutionary shocks are required only 
on the boundaries between domains in parameter space that contain only evolutionary 
shocks. 

6.2. Steepening of continuous waves 

Another argument that appears to justify the existence of intermediate shocks is based 
on the results of numerical simulations by Wu (1987), which suggest that intermediate 
shocks can be formed by nonlinear steepening of simple magnetosonic waves. Since the 
transverse component of magnetic field changes sign across an intermediate shock the 
simple wave must have the same property, which means that the transverse component 
of magnetic field must vanish somewhere within the wave. However, at this point the 
magnetosonic speed is equal to the alfven speed and it is impossible to assign a unique 
eigenvector to the simple wave. As the result, the direction of the tangential component of 
the field can rotate by an arbitrary angle at this point so that simple wave really consists 
of two distinct parts, which are disconnected as far as the direction of the magnetic field 
is concerned. This can be put in a slightly different way. Alfven waves propagating in the 
same direction as such a simple wave cannot pass throught the alfven point. During the 
steepening they will accumulate near this point giving rise to a net field rotation so that 
the discontinuity that forms has non-coplanar left and right states and can therefore not 
be a single shock. Instead, it must split into evolutionary shocks, one of which must be 
an alfven shock. Incidentally, this seems to be the only way of generating alfven shocks. 
However, in planar MHD the transition through the alfvenic point is unique and as 
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Figure 7. MHD shocks interacting with small amplitude alfven waves. The evolutionary shocks 
survive, but the non-evolutionary ones split, (a) A 1 — > 3 shock splits into a fast shock (FS), an 
alfven wave (AW) and a slow shock (SS); b) A 2 — > 4 shock splits into an alfven wave and a slow 
shock. Other small amplitude waves are also emitted. The dashed line shows the exact ideal 
solution of the Riemann problem for the initial state formed by the collision of the intermediate 
and alfven shocks. 



we have seen some of the intermediate shocks are in fact evolutionary. This is the 
explanation for the outcome of the planar simulations performed by Wu (1987). He also 
found that the results were not very different if the initial data was perturbed so that it 
was no longer exactly coplanar. However, because of the periodic boundary conditions 
used in this simulation, there was no net rotation in the perturbed problem, which makes 
it rather artificial. The reason why this perturbation did not destroy the intermediate 
shock is that these boundary conditions, together with the initial data, only allowed a 
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Figure 8. Brio & Wu problem (Brio & Wu (1988). (a) Numerical solution found using a 
Godunov type scheme. This is a proper solution of the reduced system of planar MHD but is 
inadmissible in full MHD. (b) Numerical solution found using Glimm's scheme to track alfven 
discontinuties (markers) and the exact solution involving only evolutionary shocks (lines). This 
is a proper solution for full MHD and is the only physically admissible solution for this problem. 

small value of I z per shock. It is therefore hardly surprising that an intermediate shock 
appeared since, as we have shown, these shocks can survive if I z is small enough. 

6.3. Timescale for disintegration 

Let us suppose that an intermediate shock has somehow been formed and then interacts 
with an alfven wave that rotates the magnetic field by a small angle 6<j>. It is clearly 
of some importance to know how long it takes for the shock to split. Our simulations 
show that it splits when the the value of I z associated with the shock structure becomes 
comparable with lB y , where I is the shock thickness. If the incident alfven wave has a 
small amplitude, 6(j>, then this gives us the following estimate for the disintegration time, 
t. 



I 



(6.1) 
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where we have used the alfven speed as a characteristic fluid velocity in the shock frame. 
This also tells us that the shock will only propagate for a distance & l/5<j) before it falls 
apart. From this we conclude that in all cases for which the dissipative scale is much 
smaller then the characteristic length scale of the flow, intermediate shocks can only 
appear as very short lived time-dependent phenomena. 

It is instructive to apply equation (3.1) to the interplanetary intermediate shock for 
which Chao et al. (1993) claim to have found evidence in the Voyager 1 data. In this 
case c a — 40 km s _1 and I = 5 x 10 4 km, which gives t s — 1.210 3 <5</> _1 s. The flow time 
for the solar wind at this distance (~ 9 AU) is ~ 3 x 10 7 s. It is therefore clear that Scf> 
would have to be ridiculously small for the shock to survive for a significant fraction of 
a flow time. This is most unlikely since the flow of the solar wind is sufficiently complex 
to contain plenty of alfven waves for which 5(f) ~ 1 and indeed Chao et al. find plenty of 
evidence for strong alfven waves in the data. Actually, the evidence for an intermediate 
shock is not really very convincing. The uncertainties are such that it could just as well 
be a slow shock. 

Exactly the same arguments can be applied to the magnetohydrodynamic shocks in 
the interstellar medium. Not only does the theory of collisionless shocks (see e.g. Tidman 
& Krall 1971) predict that, in these conditions, such shocks are extremely thin compared 
to the scale of the flow but there are numerous observations that confirm that this is 
indeed true (see e.g. Draine & McKee 1993). 



6.4. Convexity of MHD 

From the above discussion it is quite clear that a hyperbolic system is genuinely non- 
convex if it allows structurally stable compound waves that only contain evolutionary 
shocks. Planar MHD is therefore genuinely non-convex whereas full MHD is convex. 



6.5. N on- evolutionary shocks in numerical simulations 

The appearance of non-evolutionary shocks in numerical calculations is not something 
that is unique to MHD since it is well known that, even in gas dynamics, some numerical 
schemes can generate expansion shocks in certain circumstances. However, this phe- 
nomenon is both more subtle and more interesting in the case of MHD. The essential 
point is that, unlike gas dynamics, planar MHD is is very different from full MHD in 
the sense that there are shocks that are non-evolutionary in full MHD, but evolutionary 
in planar MHD and vice-versa. Unfortunately, this property means that the results of 
planar MHD simulations can be very misleading because, although most upwind schemes 
seem to give perfectly good solutions for planar MHD, these are of no relevance to the 
real universe with its three spatial dimensions. This is not at all unusual, indeed it may 
very well be the rule rather than the exception. For example, the properties of fluid 
turbulence are very different in two and three dimensions as are those of magnetohydro- 
dynamic dynamos. 

The other properties of non-evolutionary MHD shocks, that are not shared by gas dy- 
namical expansion shocks, are that all of them satisfy the second law of thermodynamics 
and most of them also possess a steady dissipative structure. This, together with the 
fact that the ratio of the thickness of numerical shock structures to the overall scale of 
the flow is almost always many orders of magnitude greater than in the corresponding 
physical system, means that they can persist for a significant time even in nonplanar 
problems. For example, if the piston problem discussed by Jeffrey & Taniuti (1964, p. 
256-258) is slightly modified so that it has a small transverse component of the field, 
then the evolutionary solution contains fast, slow and alfven shocks all propagating with 
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very similar speeds. In a numerical simulation this complex would remain unresolved for 
some time, during which it would be classified as a 1 — * 4 shock. 

The only truly satisfactory solution to this difficulty is to devise schemes that only 
allow evolutionary shocks. Figure 8 shows that there are schemes that will do this. 
Here we have a numerical solution to the Brio & Wu problem obtained with our MHD 
version of Glimm's scheme (Glimm 1965). This method requires a nonlinear Ricmann 
solver and we employ the one described in Falle et al. (1998), which specifically excludes 
intermediate shocks. In fact we do not use Glimm's scheme everywhere, but only to 
track the alfven shock. One can see that in this way we can avoid the appearance of 
intermediate shocks even in planar problems. Unfortunately, it is not a simple matter to 
generalise this to more than one dimension. 

The only viable option, that we can think of, is to subject all numerical calculations 
to a careful analysis using the theory described in this paper. As an example of this, it 
is instructive look at some recent calculations of steady MHD flow past a cylinder. 

6.6. 2D bow shock simulations 

Recently, De Sterck et al.(1998) have carried out numerical MHD calculations of the flow 
past an infinite, perfectly conducting cylinder. These are planar simulations and must 
therefore be interpreted in the light of the theory of planar MHD. The parameters are 
chosen in such way that the usual convex bow shock is impossible. Instead, the analysis 
given in Steinolfson & Hundhausen (1990) suggests that the shock has a dimple. They 
assumed that there is only a single shock, in which case a cosistent solution requires the 
shock type to change from 1 — > 2 to 1 — > 3 and then to 1 — ► 4 as the distance from the 
symmetry axis decreases. Although the 1 — > 4 shock is non-evolutionary even in planar 
MHD, in this case it seems that such a shock must occur on the symmetry axis for the 
same reason that it occurs when a piston moves parallel to the magnetic field. However, 
one would expect it to split into 1 — > 2 and 2 — s- 4 or 1 — > 3 and 3^4 shocks further 
away from the the axis. Indeed, De Sterck et al.(1998) find that not far from the axis 
the l—>4 shock splits and the leading shock (ED in their notation) is a 1 — > 2. At some 
distance from this branching point the other shock (EG) is identified by them as / — s- s, 
but this is unlikely to be true everywhere for such an inhomogcneous flow. One would also 
expect another branching at the point where Steinolfson & Hundhausen (1990) predict 
a transition from 1 — > 3 to 1 — > 4. The results of De Sterck et al.(1998) do, indeed, show 
this branching (DE and DG), with the trailing shock being clearly identifiable as a 2 — > 4 
shock. 



7. Conclusions 

Both our analysis and numerical results show that the evolutionary conditions for 
existence and uniqueness of discontinuous solutions of the equations of ideal MHD are 
not only compatible with the conditions for existence and uniqueness of steady dissipative 
shock structures, they are actually complementary to them. The general theory suggests 
that this will be true for all nonlinear hyperbolic systems that can arise in nature. Non- 
evolutionary shock can have a nonuniquc dissipative structure and may, perhaps, appear 
under some exeptional curcumstances as transient phenomena. However, they are not 
persistent and are bound to split when subjected to small perturbations. In the case of 
MHD, alfven waves are the most effective killers since not only our calculations but also 
those described by Wu (1988a) show that intermediate MHD shocks are destroyed by 
interactions with alfven waves. It is true that it takes a finite time for this interaction 
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to take place, but in any physical system that we know of, this time is so short that it is 
most unlikely that such shocks can be detected. 

The occurrence of intermediate MHD shocks in planar numerical simulations is con- 
sistent with the mathematical properties of planar MHD, in which 1 — > 3 and 2^4 
shocks become evolutionary but the alfven shock becomes non-evolutionary. However, 
the planar limit is a singular limit of full MHD and we suggest that planar numerical 
simulations should be avoided, especially since they are hardly any cheaper than for full 
MHD. 

Intermediate shocks may even pollute full MHD simulations because numerical shock 
structures are usually not very thin compared to the length scale of the flow. It is 
therefore essential that the results of such simulations be subjected to a careful analysis 
in order to make sure that they do not contain any intermediate shocks. If they do, then 
additional work is required to determine the extent to which they are corrupted. The 
results of our calculations with Glimm's scheme show that this problem can be eliminated 
in numerical schemes that treat shocks especially alfven shocks, as discontinuities. 
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